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A finite element method using B-splines is presented and compared with a conventional fi¬ 
nite element method of Lagrangian type. The efficiency of both methods has been investigated 
at the example of a coupled non-linear system of Dirac eigenvalue equations and inhomogeneous 
Klein-Gordon equations which describe a nuclear system in the framework of relativistic mean field 
theory. Although, FEM has been applied with great success in nuclear RMF recently, a well known 
problem is the appearance of spurious solutions in the spectra of the Dirac equation. The question, 
whether B-splines lead to a reduction of spurious solutions is analyzed. Numerical expenses, preci¬ 
sion and behavior of convergence are compared for both methods in view of their use in large scale 
computation on FEM grids with more dimensions. A B-spline version of the object oriented C-|—I- 
code for spherical nuclei has been used for this investigation. 

PROGRAM SUMMARY 


Title of program-. bspFEM.cc 

Catalogue number: . 

Program obtainable from: 

Computer for which the program is designed and others on which it has been tested: any Unix work-station. 

Operating system: Unix 

Programming language used: C-|—I- 

No. of lines in combined program and test deck: 

Keywords: B-splines, Finite Element Method, Lagrange type shape functions, relativistic mean-field theory, mean-field 
approximation, spherical nuclei, Dirac equations, Klein-Gordon equations, classes 


Nature of physical problem 

The ground-state of a spherical nucleus is described in the framework of relativistic mean field theory in coordinate 
space. The model describes a nucleus as a relativistic system of baryons and mesons. Nucleons interact in a relativistic 
covariant manner through the exchange of virtual mesons: the isoscalar scalar cr-meson, the isoscalar vector a;-meson 
and the isovector vector p-meson. The model is based on the one boson exchange description of the nucleon-nucleon 
interaction. 


Method of solution 

An atomic nucleus is described by a coupled system of partial differential equations for the nucleons (Dirac equations), and 
differential equations for the meson and photon fields (Klein-Gordon equations). Two methods are compared which allow 
a simple, self-consistent solution based on finite element analysis. Using a formulation based on weighted residuals, the 
coupled system of Dirac and Klein-Gordon equations is transformed into a generalized algebraic eigenvalue problem, and 
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systems of linear and nonlinear algebraic equations, respectively. Finite elements of arbitrary order are used on uniform 
radial mesh. B-splines are used as shape functions in the finite elements. The generalized eigenvalue problem is solved 
in narrow windows of the eigenparameter using a highly efficient bisection method for band matrices. A biconjugate 
gradient method is used for the solution of systems of linear and nonlinear algebraic equations. 

Restrictions on the complexity of the problem 

In the present version of the code we only consider nuclear systems with spherical symmetry. 

LONG WRITE-UP 


I. INTRODUCTION 

Over the last decade, the relativistic mean field theory (RMF) has been applied with great success 
to the description of low energy properties of nuclei gjlj and to the description of scattering at 
intermediate energies []. Therefore, RMF gains increasing recognition. Effective models have been 
suggested which are represented by Lagrangians containing both, nucleonic and mesonic fields 
with coupling constants that have been a(Rusted to the many body system of nuclear matter and to 
finite nuclei in the valley of /3-stability Of course, such a procedure is completely phenomeno¬ 
logical and in spirit very similar to the non-relativistic density dependent HF-models (DDHF) of 
Skyrme and Gogny |^Q. Compared to DDHF theory, the relativistic models seem to have impor¬ 
tant advantages: (i) they start on a more fundamental level, treating mesonic degrees explicitly and 
allowing a natural extension for heavy-ion reactions with higher energies, (ii) they incorporate from 
the beginning important relativistic effects, such as the existence of two types of potentials (scalar 
and vector) and the resulting strong spin-orbit term, a new saturation mechanism by the relativistic 
quenching of the attractive scalar field and the existence of anti-particle solutions, (iii) Hnally they 
are in many respects easier to handle than non-relativistic DDHF calculations. 

Since the discovery of the halo phenomenon in light drip-line nuclei ^ the study of the structure 
of exotic nuclei has become a very exciting topic. Experiments with radioactive beams provide a lot 
of new data over entirely new ("exotic”) regions of the chart of nuclides. On the theoretical side, 
presently existing models of the nucleus, relativistic ones as well as non-relativistic ones, have to be 
tested in these new regions in comparison with experiment. Improvements and extensions of the 
models become necessary. 

Recent investigations ]9|jic[| have shown that coupling to the particle continuum and large exten¬ 
sions in coordinate space have to be taken into account in order to describe phenomena of exotic 
nuclear structure. The underlying equations of all nuclear models have therefore to be solved on 
discretizations in coordinate space. In contrast to "conventional” methods, based on expansions of 
the solution in basis functions with spherical or axial symmetry, sophisticated techniques have to 
be applied in order to solve the mean field equations in coordinate space. 

With the non-relativistic HF-models extensive nuclear structure calculations have been performed 
based on the imaginary time method [^. This very efficient method, however, is restricted to 
the non-relativistic cases where the single particle spectrum is limited from below. In relativistic 
model calculations, the imaginary time method would not converge due to mixing with negative 
energy states. Therefore, we plan a different approach with Krylov-subspace based methods 
(for solutions on 2D and 3D meshes in coordinate space) and with the bisection method (ID spherical 
case 1^]). In contrast to the imaginary time method, the required single particle or quasi particle 
eigenstates have to be calculated in each step of a self-consistent iteration. At hrst sight, it seems, 
that this approach is intractable since coordinate space discretizations on 2D or 3D finite element 
meshes lead to eigenvalue problems of large dimensions. With the block Lanczos method however, 
the calculation of eigenvalues and corresponding eigenvectors can be restricted to a small number 
which is required in the region of bound nucleons. In combination with the selfconsistent iteration 
method which is applied to the whole problem, the number of internal block Lanczos iterations can be 
reduced to corrections of the vectors which come from the previous iterations step of the selfconsistent 
loop. In references [ p'slJl^ the solution of the spherical RMF equations and the spherical RHB 
equations with the hnite element method in coordinate space has been demonstrated. In these 
investigations, I have observed that spurious solutions appear in the spectrum of eigenvalues of the 
Dirac operator of the RHB equations when they are discretized with hnite elements of the Lagrangian 
type. Since the numerical cost to calculate eigensolutions on ID-meshes is relatively small, it was not 
important to avoid spurious solutions a priori and therefore they have been eliminated by comparison 


2 


of the number of nodes. In the 2D and 3D cases, however, it is important to reduce the size of the 
stiffness matrices to a minimum. This can be achieved by using shape functions with extremely 
good properties of interpolation, allowing wider meshes in coordinate space. Since B-splines are 
smooth, one would expect that they have the desired properties. 

The major goal of the present paper is to give an answer to the question whether B-splines can 
improve the numerics in comparison to the often used shape functions of Lagrangian type. At the 
present state, our study is restricted to the solution of relativistic mean field equations. The results 
of our investigation are important with respect to large scale computations on finite element meshes 
of two and three dimensions. Such calculations are required in the relativistic mean field description 
of deformed exotic nuclei at low energies. I have worked out a B-spline version of the computer code 
which is published in |j^ and compare the results obtained with both codes for spherical nuclei. 


II. THE RELATIVISTIC MEAN FIELD EQUATIONS 


The relativistic mean field model describes the nucleus as a system of nucleons which interact 
through the exchange of virtual mesons: the isoscalar scalar (T-meson, the isoscalar vector Ci;-meson 
and the isovector vector p-meson. The model is based on the one boson exchange description of the 
nucleon-nucleon interaction. The effective Lagrangian density is ^ 


jC = Ip {i'y ■ d - 

+ liOaf 

-gaipcnp 


m) ^p 

- U{a) - -f -f im^p" 

- Qui'ipy ■ i^'tp — Qp'ip'y ■ pr'tp — eipj ■ ’ 



Vectors in isospin space are denoted by arrows. The Dirac spinor tp denotes the nucleon with mass 
m. rua, rriui, and nip are the masses of the a-meson, the u-meson, and the p-meson. ga, Qui, and gp 
are the corresponding coupling constants for the mesons to the nucleon, e^/dvr = 1/137.036. 

Since the relativistic mean field model has been described in a large number of articles, I omit a 
long discussion of the above given Lagrangian and the derivation of the RMF equations. Instead, 
I refer to section 2 of reference and to section 2 of reference [Q. In these references, the 
development preceding to the investigations of the present paper is described in details. The main 
interest of the work presented below is focused on numerical aspects and performance of two FEM 
techniques in the solution of the RMF equations for spherical nuclei. In the following, I briefly list 
the static RMF equations for the spherical symmetric case. 

Introducing spherical polar coordinates (r, 6, (p), the Dirac equation reduces to a set of two coupled 
ordinary differential equations for the amplitudes g{r) and f{r) for proton and neutron states 


+ (m-l- S'(r) - V{r)) f{r) = -e fi{r), 
fir) + {m + S{r) -|- V'(r)) 5 i(r) = +eg{r), 


where the quantum number k = ±1, ±2, ±3,.... The scalar potential S{r) and the vector potential 
V (r) are composed of boson field amplitudes and coupling constants where 

5(r) = srao-(r). 


and 

V{r) = g^uJ°{r)+gpT 3 p° 3 {r)+e ^^ A°(r). 

The symbols g^, gui, gp, and e denote the coupling constants of the cr-field, the cu-field, the p-field 
and the A-field, coupled to the nucleons. The meson fields a{r), uj^{r), Psir) and the photon field 
A°(r) are solutions of the inhomogeneous Klein-Gordon equations 

(^-dr - ^dr + + ml'j a{r) = -g^ Ps(r) - 52 cr^(r) - ps cr®(r) 

(^-dr - ^dr -I- -I- mlj 01 °(r) = g,^ p^{r) 


( 1 ) 


( 2 ) 


(3) 

(4) 


(5) 

( 6 ) 
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(^-dl - p°(r) = Qp psir) 

(^-dr - ^dr + A°(r) = epeni(r) 

where the sources of the fields are the scalar density Ps{r), the isoscalar baryon density pv{r), the 
isovector baryon density P3(r) and the electromagnetic charge density. They are composed of the 
nucleon wave functions in a bilinear way as 

P^{r) = ^ n^,n^^^{gK,r,{rf - fn,n{rf) 

K, n 

pv = ^ + fn.nirf'j 

K,n 

K,n 

Pem = ^ 

K,n 

where the quantities nn,n are occupation numbers of the energy levels (indices Ac,n). For the simple 
Hartree model without pairing, = 1 for occupied levels and = 0 for unoccupied levels. 
The index n denotes the principal quantum number (n = 0,1,2,...) and counts the eigensolutions 
of equation from small to large energies n- The nucleon numbers filling an orbital (K,n) are 
taken into account by the factors 2\k\ in Eqs. Since the densities (|)-® do not depend 

on the aMular coordinates 9 and (j), no terms higher than of monopole order show up at the r.h.s. 
of Eqs. (^ - (H)- Consequently, the solution of the these equations has to be restricted to Z = 0 in 
the description of spherical nuclei. For physical reasons, the nonlinear self-coupling of the n-field 
has to be taken into account. It is described by the two terms —g 2 and —gs a-^{r) at the r.h.s. 

of Eq. (^. Without these terms, the RMF-model could not explain the compressibilities in finite 
nuclei and nuclear matter as well as the surface properties of finite nuclei. 


(7) 

( 8 ) 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 


III. B-SPLINE AND LANGRANGIAN TYPE FINITE ELEMENTS 


The most widely used finite element type in many applications is the Lagrange type element. 
Lagrangian shape functions allow the simplest representation compared to other types of shape 
functions. For any finite element order n they have the following expression in reference element 
representation (see Figs. 2a, 2c, 2e) 


n 


Kip )=n 


z=o 


jnp - 1) 

(i-O 


where the coordinate p is restricted to the interval [O, l]. From Eq. it becomes obvious that 
Lagrange type finite elements are easy to handle. Generally, for any conventional finite element 
type, the shape functions (nodal basis) have the property 


KWM = V 


in the one dimensional case and 


ipg') ^qq' 5 P — (Pli Pm) 

in the M-dimensional case where Pgi denotes a grid point of M-dimensional finite elements. These 
shape functions form a so called nodal basis. Shape functions with property ( p^ can be constructed 
for finite elements of any geometrical form as for example triangular elements or quadratic elements 
in two dimensions and tetraedric or cubic elements in three dimensions. Also the location of the 


(13) 


(14) 


(15) 
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mesh points which belong to a finite element can be distributed in almost arbitrary manner over the 
domain of the element. In most cases, however, M-cube elements (intervals, squares, cubes, etc.) 
with a uniform distribution of the nodes are sufficient and allow extremely efficient calculations of 
the stiffness matrices for a given boundary value problem. The shape functions of such elements are 
represented as products of Lagrange polynomials ( p^ 

M 

(f'l’ ■ • ■, pm) = iVq"; (pi), qi,:;ni, 

i=\ 

where rii denotes the order of the element in the direction of dimension i and {q) = {qi, qM) 
forms the index tuple of the nodes. 

The construction ( [l^ of Lagrange type M-cube shape functions shows that 1-dimensional La¬ 
grange type finite elements allow the most simple generalization to M-cube meshes. A great advan¬ 
tage of such shape functions can also be seen from a technical view point. Implementations of the 
general M-dimensional case in object oriented programming styles become simple. The amount of 
data required by an object which represents a Lagrangian M-cube finite element as reference element 
is almost the same as in the 1-dimensional case if the orders ni,...,nM are equal. In this special 
case the data representing all shape functions of a M-cube element comprise 2 ■ rii • nf floating point 
numbers where I denote by nf the number of Gauss points on a Gaussian mesh in dimension i. In 
the more general case where m yf nj for i j and nf yf nf for i j the amount of float point 
numbers required to represent all shape functions is 

M 

2 ^ m nf 

i=l 

which is still small compared to the number of values 


M 

2 P[ ni nf 

i=l 


required for shape functions of arbitrary type and arbitrary distribution of the nodes over the 
element. 

The numerical cost for the integration of matrix elements reduces dramatically in cases where 
operators split up into products of operators each depending on a complementary subset of the 
coordinates pi,...,pM- In the most ideal case an operator factorizes completely leading with Eq. 
(0) to a complete factorization of matrix elements. 


IVl IVl 

..«*')) = 


Oi{pi 


Knp^))- 


i=l i=l 

This becomes obvious when I rewrite Eq. ( p^ in terms of a numerical Gauss integration 

N M M N 




i=l 


i=l L=l 


where, assuming that N is the number of Gauss points in each coordinate, the number of floating 
point operations on the left hand side is greater than 4N^ while on the r.h.s. it is smaller than 


3MN. 


( 16 ) 


(17) 


(18) 


(19) 


( 20 ) 
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Fig. 1: B-spline functions of different orders increasing from 1 to 5. 


These advantages of Lagrangian type finite elements gave us a reason to apply them in several 
previous studies and calculations (see references In these references, it has been shown 

that Lagrangian finite elements provide an excellent tool for solving the equations of the RMF 
model in self-consistent iterations in coordinate space. In this manuscript, I present a new finite 
element technique using B-splines as shape functions and compare this method with the FEM based 
on Lagrangian shape functions. B-splines have a compact support and are defined as polynomials 
piecewise on intervals which are bounded by neighbored mesh points. The basic criterium in the 
construction of these basis functions is optimized smoothness over the whole support. This property 
is guaranteed if all derivatives up to the order n — 1 obey the conditions of continuity in all matching 
points of the mesh. The order n of a given B-spline corresponds to the degree of the polynomials 
by which it is composed. In Fig. 1, examples are shown for B-splines of order one to five, n + 2 
mesh points are required to construct a B-spline of the order n. In contrast to Lagrangian shape 
functions, B-splines of any order do not change sign. A common property of both types of shape 
functions is expressed by 

p p 


where p denotes the mesh point index. The fact, that B-splines of any order satisfy all these 
conditions makes it impossible to find an expression in closed form in the sense of Eq. ( |l3| ) . Rather, 
they are generated by the following brief algorithm. 


start: 




{{Xi+l - Xi)) ^ Xi < X < Xi+l 

0 X < Xi, X > Xi+i 


Bi,k{x) = 


{X - Xi) 
{Xi+k - Xi) 


Bi^k-i{x) + 


{Xj+k - x) 
{Xi+k - Xi) 


Bi+i^k-i{x) 


( 21 ) 


( 22 ) 
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I define a B-spline finite element as a region which is bounded by two neighboring mesh points in the 
one-dimensional case or as a M-cube where the 2^ corners are identical with the 2^ mesh points 
of a cubic grid which are closest to the center of the cube. Obviously, this definition is restricted 
to cubic grids but it will turn out to be extremely efficient in all cases where cubic finite element 
discretizations can be applied. 

The figures 2b, 2d, and 2f display 3'^'^ order, 4**“ order and S*** order B-spline finite elements in one 
dimension according to our definition. The figures show, that all parts of a B-spline are contained 
in each element. 

Fig. 2a Fig. 2b 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 " Vo 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



Fig. 2 : Finite elements in reference element representation. In the figures (a),(c),(e), Lagrange 
type elements with corresponding shape functions are shown for 3'^'^ order, 4*^ order and 5*^ order. 
In comparison, B-spIine elements of 3“^*^ order, 4*^ order and S*** order are displayed in the figures 

(b), (d), (f). 
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A B-spline of order n is determined by n + 2 mesh-points through the algorithm and the 
support is composed by n -|- 1 intervals. However, the pieces which are contained in a single element 
belong to different B-splines each located at another mesh point and determined by a different set 
of mesh points. They are the overlap sections of all splines which are not zero in the considered 
element. This fact has to be taken into account especially if nonuniform meshes are used. As a 
consequence, the amplitudes of the shape functions in such a finite element depend on the position 
of all mesh points which are in the support of all contributing B-splines. Most of these mesh points 
are outside of the element and belong to other elements. Therefore, the degree of ’’interaction” 
between neighbored elements is maximized and much higher than for elements of Lagrange type 
where only next neighbor elements ’’interact” through their outermost contact grid points. This 
leads to a strong reduction of the total number of required mesh points as I will demonstrate in 
section 5. 

A disadvantage, however, appears when non-uniform meshes are used. As it has been discussed 
in reference |^], Lagrangian elements can be mapped on a reference element through linear affine 
transformations. This is even possible for non-uniform meshes. Therefore, Lagrangian shape func¬ 
tions need to be evaluated only once in the reference element and amplitudes and abscissas can 
be accessed by means of a pointer to that reference element. These advantages are also valid for 
B-splines as long as uniform meshes are used. In the case of non-uniform meshes the algorithm 
(^ ) has to be evaluated for each argument taken on the global region. This leads to a reduction 
of storage requirement but increases the numerical cost by a factor which corresponds to the num¬ 
ber of floating point operations which are necessary to carry out the scheme (M). Consequently, 
the numerical effort in the calculation of the stiffness matrices of a given problem increases by the 
same factor. On the other hand, the number of algebraic equations resulting from a finite element 
discretization is usually large and the numerical cost to solve these equations increases faster with 
the number of equations (~ number of mesh points) than the numerical cost in the calculation of 
the stiffness matrices with the number of mesh points. This trend is even enhanced for increasing 
dimension M of the descretization where the condition of the stiffness matrices becomes worse. 
B-splines finite elements might therefore be superior in multidimensional FEM discretizations as 
compared to Lagrangian finite elements. 


IV. THE FEM DISCRETIZATION 


A basic principle of FEM is the approximation of the solution for a given problem in a space of 
shape functions which have compact support and existing continuous weak derivative of maximum 
degree m. Together with a p-norm which is for all those functions defined as 

^ P 1/p 

M\m,p ■- , 

a=0''^ 

the above given space is a Banach space. According to the norm it is usually called Sobolev space 
W^{Q,) where fl may be any compact domain of the coordinate space. Since W^{Q) is complete, the 
solutions of any partial differential equation of an order not higher than m can be approximated to 
arbitrary precision in W^iQ.) on the whole domain fl. This property plays an important role for the 
solution of differential equations with computational methods because finite element discretizations 
of the domain fl correspond to subsets of linear independent functions of W^{Q.) and because the 
representation of functions of on the computer is simple. In FEM, 12 is subdivided into 

a large number of small sub-domains which are called finite elements. Each element is support of 
a certain number of shape functions which is equivalent to the number of constraints set on the 
element. These functions span up a finite element space. The corners of the elements are located on 
a finite element mesh. However, additional mesh points can exist in the interior or on the surface 
of each element and additional constraints as derivatives of any order can be applied. In such cases 
the order of the element is higher than first order. 

In the following, I discuss the finite element discretization of the Eqs. (P) and Eqs (^) to for 
both types Lagrangian and B-spline elements. In the present application a nodal basis {Ap(r)} is 
used in the case of Lagrange elements and a non-nodal basis {i3p(r)} is used in the case of B-spline 
elements. Examples for discretizations with elements of order are displayed for both types in Fig. 
3a and 3b. Each Lagrangian element in Fig. 3a has two boundary nodes and two additional nodes 
in the interior whereas the B-spline elements in Fig. 3b are free of interior nodes. In the B-spline 


(23) 
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method additional nodes are required outside of the region of integration in order to generate the 
shape functions which have non-zero overlap with the inner region. I use the notation Np{p) and 
Bp{p) {0 < q < n) for shape functions in reference element representation, and Np{r) or Bp{r) 
(1 < P < for shape functions Np and Bp on the global mesh in the r-coordinate space. Using 

the standard representation for the Pauli matrices, the Dirac equation (|^ is written in matrix form 


{dr + r 


CT 3 CT 1 — Kr ^ ■ CTi -I- (m -I- S{r)) ■ az + V{r) ■ I 2 


'l>(r) = £ I 2 <l>(r). 


( 24 ) 


Fig. 3a 


Fig. 3b 




Fig. 3: Global discretizations with 3’^'^ order finite elements (a) of Lagrangian type and (b) 
of B-spline type. In both examples a total number of 16 mesh points is used in the region between 
r = 0 fm and r = 10 fm. 


For the nucleon spinor I use the FEM ansatz 


with B-spline functions Bp{r) and node variables Up ~ {Up\Up{f)Y'. In a Lagrangian FEM 
ansatz, the shape functions Bp{r) defined in (EM are replaced by the shape functions Np{r). In the 
weak formulation of the eigenvalue problem the weighted residual (see ||^) leads to algebraic 
equations of the form 

^2(wpi(r) {dr + r~^) ■ asai - Kr~^ ■ a-i + {m + S{r)) ■ as+ V{r) ■ I 2 Bp{r)'jup = 


e(wp>{r) 


Bp{r))u 


with weighting functions Wpi{r). The weighting functions are chosen 

Wpi{r) = ~ f ^ ^ Bp{r) 

\ V^max/ / 

in the case of B-spline elements or 

Wpi{r) = ~ f ^ ^ ro/1 Np{r) 

\ \'f 'max / / 

when Lagrangian shape functions are used in ansatz (|25|). 


—K — 1 K < 0; 
K K > 0; 


(25) 


(26) 


(27) 


(28) 


(29) 
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. _ J —K K < 0; 

^ — 1 K>0; 

The factor corresponds to the Jacobi determinant of the transformation into spherical coordinates 
and compensates singularities in the operator. The factor r*® or r^f respectively includes boundary 
conditions at r = 0 properly for upper {g{r)) and lower components (/(r)) of the spinor <l>(r). The 
factor (1 — (r/rmax)^) includes boundary conditions g{rma.x) — 0 at r = fmax in matrix elements 
which are multiplied with g-components in <E>(r) whereas this factor is replaced by 1 when a matrix 
element is multiplied with node variables of the /(r)-component. 


( 30 ) 



Fig. 4b 



Fig. 4d 
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Fig. 4f 




Fig. 4 : Occiipation patterns of stiffness matrices which result from the finite element dis¬ 

cretization of Eq. (0). The figures (a), (c), (e) display matrices which correspond to 3'^'^ order, 4*^ 
order and S**' order Lagrange type FEM. Figues (b), (d), (f) display patterns which result from 
order, 4*^ order and 5*^ order B-spline FEM discretizations. The numbers represent counter indices 
used in a vector storage technique. 


The boundary conditions at r 
following way. 


Ofm depend on the quantum number k and are defined in the 


g{r = 0) = 0 


0 

II 

0 

for AC = — 1 

(31) 

0 

II 

0 

for AC = -Fl 

(32) 

and f{r = 0) = 0 

for 1 ac| > 1 . 

(33) 


The system of algebraic equations ( p6| ) forms a generalized eigenvalue problem of the form Au = 
e N u with stiffness matrices A and N can be analyzed from the resulting matrix equation 


® CTscTi + S2 (Si uacri — KS2 ® cri -)- mSs (S 0-3 + S4 ^ as + S5 (S I2 


S 3 (S I2 


where u is a vector with components {u\^\uY\ . The symbols Si to S 3 denote the 

stiffness matrices of the operators on the l.h.s. of Eq. wM, 


51 = (wpi{r) dr Bp{r)'j, 

5 2 = (wpi{r) r~^ Bp{r)y 

5 3 = (wpi{r) 


Bpir)), 


S{r) 


Bp{r) 


5 4 = (wpi{r) 

5 5 = (wpi{r) V{r) Bp{r)y 


In Figs. 4a-f, occupation patterns of the stiffness matrices A are displayed for Lagrangian and B- 
spline finite element discretizations. The matrices in Figs. 4a, 4c, 4e result from the Lagrange FEM 
with order, 4*^ order and S*** order finite elements. For comparison, I show the corresponding 
stiffness matrices of the B-spline FEM in Figs. 4b, 4d, 4f. The sub-block structure of 2 x 2-blocks in 
all matrices results from the fact that Eq. (0) is a system of two coupled equations. The number of 
occupied 2 x 2-blocks for a given order in the Lagrange FEM is nB ■ + 2u] -F 1 while in 

the B-spline method the occupation increases to nB ■ |^2(n”‘^)^ -F n] -F 1. nB denotes here for both 
cases the number of finite elements used in the Lagrange method and is different from the number 
of elements which is used in the B-spline FEM. 


(34) 


(35) 
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The FEM discretization of the Klein-Gordon equations (^) to (^) is described in the appendix. 
Finally, the coupled system of differential equations (|^, (^) to (^) is replaced by a system of linear 


algebraic equations 

, A^) u = £ N u (36) 

for the node variables of nucleon spinors, and 

Ba{So\d) 3 = (37) 
= (38) 

Bp p° = (39) 

Ba A° = (40) 


for the node variables Up, ujp, pp, and Ap of the meson fields o'(r), t<J°(r), p°{r), and photon field 
A°(r). The occupation patterns of the matrices B^, B^j, Bp, and Bao for various shape functions 
are very similar to those of the matrix A (Figs. 4a-d). The main difference is that 2 x 2-blocks have 
to be replaced by single matrix elements. 


V. ANALYSIS OF SPURIOUSITY 

The appearance of spurious solutions in applications of FEM is a well known problem in general. 

First applications of the finite element method in relativistic nuclear physics [[13 have shown that 
spurious solutions appear in the spectra of the Dirac equation. Linear finite elements have been 
used to calculate solutions of the relativistic nuclear slab model. Comparisons of the solutions with 
solutions that have been obtained with other numerical techniques (shooting method) have shown 
that FEM reproduces the physical spectra very well and that spurious solutions have no influence. In 
a further step Q] Lagrangian finite elements of 1. to 4*^ order have been used in the self-consistent 
solution of the RME equations of sperical nuclei. In these calculations, it has been shown (up to 
4*^ order) that the density of spurious solutions in the spectra decreases for increasing order of the 
elements. 

In this sections, I present a systematic study of the spurious spectra which appear in the spherical 
symmetric case. In the initial step of a self-consistent ground state calculation of IfPb, Woods- 
Saxon potentials 

^(r) = S'(0)(^1-bexp(^^-^))^ , (41) 

Y(r) = Y(0)(l+exp(^l^))"\ (42) 

are used for the scalar potential S{r) and for the vector potential V{r). For the values 

of these potentials at r = Ofm are chosen S'(O) = —395MeV and Y(0) = 320MeV, respectively. 
a = 0.5 fm and = 9.0 fm. The calculation is performed on a uniform radial mesh extending 

from rmin = Ofm to rmax = 20fm. A smaller value rmax = 12fm would be sufficient for 
to obtain good approximations of the bound single particle states. For a good resolution of the 
continuum, however, a large extension of the mesh in coordinate space is necessary. An extremely 
high number of 200 mesh points has been used in the calculation of the sprectra shown in Fig. 5a 
to Fig. 5f. The reason for that choice will become clear from the subsequent discussion of Fig. 

6 . For a nucleon mass of 939 MeV (parameter set NL3), the Dirac gap extends from —939 MeV 
to +939 MeV. Bound solutions are expected to have energies which are located in the Dirac gap. 

In the following calculations an energy window ranging from —1300 MeV to +1300 MeV has been 
chosen which covers parts of the lower and upper continuum as well. 

The results which are presented in the subsequent discussion correspond to the first iteration step 
and a value k = —1 (s-waves). Spurious spectra of similar eigenvalue distributions are obtained for 
all other K-values (k = +1, ±2, ±3,...). Disregarding the fact that the eigensolutions change while 
they converge, very similar results are found in all iteration steps of the selfconsistent iteration. 

One of the most interesting questions to be answered in the present paper is, whether the ap¬ 
pearance of spurious solutions can be avoided using B-spline finite elements instead of Lagrangian 
elements. Since both methods are identical in the case of F* order, spurious solutions appear also 
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in the B-spline FEM. However, from that one can not conclude that spurious solutions appear in 
FEM discretizations with B-splines of higher order. The following six figures Fig. 5a to Fig. 5f 
display energy spectra of Eq. (|^ which have been calculated for many different orders with both 
methods, the Lagrange FEM and the B-spline FEM. In Fig. 5a and Fig. 5b the positive and neg¬ 
ative spectra are shown for order to d**' order finite elements. The white circles correspond to 
physical eigenvalues which have been calculated with Lagrange type elements. They are located at 
the same energies as the white triangles which correspond to physical eigenvalues obtained with the 
B-spline FEM. The figures show that the physical spectra are independent on the order and on the 
used method. However, the number of black filled circles and filled triangles in the spectra decreases 
for increasing order of the used shape functions. The filled symbols indicate eigenvalues of spurious 
solutions. It turns out that the distributions of spurious eigenvalues over the entire energy range 
are identical for both methods and in all orders. In comparison to the Lagrange FEM, the B-spline 
method does obviously not reduce the number of spurious states as long as the order is the same 
used in both methods. However, for both methods, the density of spurious solutions in the spectra 
can be strongly reduced by increasing the order of the hnite elements. Particularly, from Fig. 5c 
to Fig. 5f, on can see, that the spurious eigenvalues drift away from the Dirac gap when the order 
is increased. Consequently, for any energy window there exists an order which is high enough so 
that no spurious solutions appear in the window. An exception forms the region between the lowest 
positive physical eigenvalue and the highest negative physical eigenvalue. For all element orders 
with both methods, no spurious solution has been found in that region. 

Fig. 5 a Fig. 9 b 




Fig. 5 c Fig. 5cl 
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Fig. 5 e 


Fig. 5 f 




Fig. 5: Energy eigenvalues of the Dirac equation (^) for the case k = —1. The spectra 

are compared for the Lagrange FEM (circles) and the B-spline FEM (triangles). The used finite 
element orders are indicated in the figures. Filled symbols correspond to spurious eigenvalues. All 
eigenvalues which appear in the energy window [—ISOOMeV, +1300MeV] are displayed for I*** order 
to 12*^ order elements. 

It is also important to analyze the dependence of the spurious spectrum on the number of mesh 
points. In Fig. 6, the number of spurious solutions is displayed for different orders as a function of 
the number of mesh points in a constant radial box (rmin = Ofm, Tmax = 20 fm). The results show 
that the number of spurious solutions is independent on the number of mesh points if this number 
is sufficiently large. This is true for all orders of finite elements. The solid lines in Fig. 6 show 
the results which have been obtained for the above defined Woods-Saxon potentials. For all finite 
element orders the number of spurious states in the above defined energy window increases at low 
mesh point numbers and decreases monotonically at high mesh point numbers. At large numbers 
of mesh points (’’asymptotic region”) the number of spurious solutions is constant for all element 
orders. This has been tested up to the very large number of 600 mesh points but is not shown in 
the figure. For the calculation of the spectra shown in Fig. 5a to 5b, I have used a number of mesh 
points (200) which is in that asymptotic region to make sure that they (in particular the spurious 
spectra) are independent on the number of mesh points. 

An explanation for the curves in Fig. 6 is found with the concept of Sobolev space. In reference 
0 it has been outlined that a Sobolev space is a completion of the test function space 

C^(fl) with respect to the Sobolev norm || • \\m,p defined in Eq. (|23|). Thus, C W^iQ.) for 

all integer numbers m > 0. All spaces where m > 0, are subspace of the largest Sobolev 

space IFp(fl) and 

W™+^(D) C IFp"*(D) for all m > 0. 

The shape functions of m*** order finite elements are element of W^{Q,) but the shape functions of 
any lower order finite elements are not in W^(D). In finite element discretizations of low order m 
is small and one works in a correspondingly large space W^iyi). The weak form, of a differential 
equation, expressed in terms of the weighted residual, allows more solutions than the solutions of 
the original problem. All solutions which are found for a certain FEM order m are element of 

For increasing order m the space W^iQ) shrinks and the number of spurious solutions in 
the weak form is reduced while all physical solutions are maintained. This is seen from Fig. 5a to 
Fig. 5f for one large number of mesh points. It explains in general the reduction of the number 
of unphysical solutions in Fig. 6 for all mesh point numbers when the order of the FEM-ansatz is 
increased. For a uniform finite element mesh with constant width h, the m*’*' order shape functions 
of the whole FEM discretization span up a space Starting from an initial discretization 

where ho is large, a sequence of spaces SJ^. (fl) i = 0, 1, 2,... is generated when the mesh is refined for 
increasing index i, where hi+i < hi. The direct sum of the spaces S^(fl) converges against W^iO) 

OO 

and thus ^ (fl) = Different spaces (fl) and (fl) where i ^ j can have non-trivial 

i —0 

intersection. There are even cases where C Sf[](Q) when hi < hj. In the example of the 

two spaces Sl^{Q) and it is obvious that Sl^{Q) C since each linear shape function 


(43) 
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which is basis function in Sl{Q,) can be represented as linear combination of shape functions (basis 
functions) of The strong increase of the graphs in Fig. 6 at small mesh point numbers, 

where h is large, is explained by the fact that the spaces SJ^ (fl) become large for decreasing h. The 
number of spurious solutions which appear in increases simultaneously. However, there is a 

second effect which is superposed to this first one. 



number of mesh points 


Fig. 6: Dependencies of the number of spurious solutions on the number of mesh points 

are shown for T** order to 6*^ order finite elements. A constant mesh size ranging from Ofm to 
20 fm has been used in the calculations. The solid (and dot-dashed) lines show results obtained 
for Woods-Saxon potentials V{r) and S{r). The dashed lines show corresponding results for zero 
potential. 

Spurious solutions which have a very high number of oscillations can only be completely resolved in 
spaces S'))*(D) where h is small. However, spurious states with high frequency can appear in subspaces 
S™(D) where h = v ■ h {u = 1,2,3 ,...) and where only a fraction of the oscillations is resolved. Since 
the corresponding kinetic energy which contributes to the total energy is small, these solutions 
appear in the above chosen energy window. If the number of mesh points is increased, additional 
oscillations are resolved and the kinetic energy increases correspondingly. The corresponding total 
energy appears no longer in the energy window. In the negative energy range such solutions are 
shifted further into the negative continuum. 

In Fig. 7a and Fig. 7b, spurious energy spectra are displayed for many different mesh point 
numbers. First order finite elements have been used. The mesh point numbers have been chosen 
around the maximum of the solid curve in Fig. 6 which corresponds to first order. The solid lines 
connect spurious eigenvalues which appear at constant mesh point number which is indicated by 
the numbers atop of each line. 
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Fig. 7a 


Fig. 7b 




index index 


Fig. 7: Spurious spectra of the first order finite element discretization are displayed for different 
mesh point numbers. The solid lines connect eigenvalues which belong to the same discretization. 
The number of mesh points is indicated above each line. 


VI. NUMERICAL PRECISION 

In the following, I present an analyse of the numerical precision of both methods and compare the 
results. The quality in the approximation of the exact solutions of (^) depends essentially of the 
order of the FEM-ansatz and on the number of mesh points used in a given domain 12 = [rniin, f’max] ■ 
In the subsequent tables, neutron single particle eigenvalues are listed systematically for increasing 
number of mesh points and increasing order of used finite elements. The eigenvalues correspond 
to solutions which have been obtained in the initial step of a selfconsistent calculation for “^^Ca. 
In the first iteration step, Woods-Saxon potentials of the form (§) and (^) have been used with 
parameters 5(0) = —395 MeV, V (0) = 320 MeV, a = 0.5 fm and rs = 6.0 fm. The mesh size has been 
kept fixed with boundaries at rmin = Ofm and r^ax = 10 fm. In Table la neutron single particle 
energies which have been calculated with linear finite elements are shown for the initial Woods- 
Saxon potential. For increasing number of mesh points (see left column), the number of unchanged 
decimal places reaches 8 at 200 mesh points. A comparison with the last row of Table lb shows 
that for linear elements the last digit (decimal place 10) has not stabilized at the extremely large 
mesh point number 600. Table lb displays results that have been calculated with finite elements of 
order. Between 109 and 121 mesh points (36-40 elements), the results have stabilized in all 10 
digits. In Table Ic, I show the corresponding results which have been obtained with finite elements 
of 4*^ order. To demonstrate the enormous improvement in the precision, the results are displayed 
up to 12 digits. Between 81 and 93 mesh points the 8*^ digit becomes stable and between 125 and 
145 mesh points the precision achieves 12 digits. 


1. Order 


'f^nod 

Eis1/2 

Elp3/2 

Flpl/2 

E]d5/2 

E2s1/2 

Eld3/2 

10 

65.88844974 

57.33594611 

56.42536540 

47.45864481 

42.33620819 

45.34520938 

20 

65.88133752 

57.31827220 

56.39907816 

47.42408909 

42.19281200 

45.28981141 

30 

65.88108502 

57.31763456 

56.39825020 

47.42279518 

42.18737854 

45.28804244 

40 

65.88104432 

57.31753150 

56.39811692 

47.42258515 

42.18648990 

45.28775645 

50 

65.88103335 

57.31750369 

56.39808101 

47.42252842 

42.18624969 

45.28767927 

60 

65.88102944 

57.31749377 

56.39806819 

47.42250815 

42.18616380 

45.28765171 

70 

65.88102777 

57.31748952 

56.39806272 

47.42249948 

42.18612708 

45.28763993 

80 

65.88102696 

57.31748747 

56.39806006 

47.42249529 

42.18610928 

45.28763422 

90 

65.88102653 

57.31748638 

56.39805866 

47.42249306 

42.18609982 

45.28763119 
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100 

65.88102628 

57.31748575 

56.39805785 

47.42249178 

42.18609441 

45.28762946 

150 

65.88102591 

57.31748480 

56.39805662 

47.42248983 

42.18608615 

45.28762681 

200 

65.88102584 

57.31748464 

56.39805641 

47.42248951 

42.18608476 

45.28762637 

300 

65.88102582 

57.31748458 

56.39805634 

47.42248938 

42.18608424 

45.28762620 

600 

65.88102582 

57.31748457 

56.39805632 

47.42248936 

42.18608412 

45.28762616 


Table la: Neutron single particle energies in units of MeV which have been calculated with linear 
finite elements of Lagrange type. The number of used mesh points has been gradually increased as 
shown in the left column. 


3. Order 


T^nod 

Eis1/2 

-Elp3/2 

Elpl/2 

Eld5/2 

E2al/2 

Eld3/2 

10 

65.87761838 

57.30204792 

56.38118314 

47.38090127 

42.16223936 

45.24653386 

19 

65.88080395 

57.31691348 

56.39762680 

47.42131836 

42.18517027 

45.28703195 

25 

65.88103067 

57.31748586 

56.39806263 

47.42248096 

42.18636302 

45.28761163 

31 

65.88103176 

57.31749678 

56.39809005 

47.42250968 

42.18614523 

45.28768121 

40 

65.88102636 

57.31748582 

56.39805841 

47.42249169 

42.18609035 

45.28762987 

46 

65.88102599 

57.31748496 

56.39805720 

47.42249002 

42.18608570 

45.28762767 

52 

65.88102587 

57.31748468 

56.39805656 

47.42248956 

42.18608460 

45.28762657 

61 

65.88102583 

57.31748460 

56.39805637 

47.42248941 

42.18608424 

45.28762624 

70 

65.88102582 

57.31748458 

56.39805633 

47.42248937 

42.18608415 

45.28762619 

79 

65.88102582 

57.31748457 

56.39805632 

47.42248936 

42.18608413 

45.28762617 

91 

65.88102582 

57.31748457 

56.39805632 

47.42248936 

42.18608412 

45.28762616 

100 

65.88102582 

57.31748457 

56.39805632 

47.42248936 

42.18608412 

45.28762616 

109 

65.88102582 

57.31748457 

56.39805632 

47.42248935 

42.18608412 

45.28762616 

121 

65.88102581 

57.31748457 

56.39805632 

47.42248935 

42.18608412 

45.28762616 


Table lb: 


Same as Table la but for Lagrange FEM of 3’^'^ 


order. 
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4. Order 


T^nod 

Eis1/2 

-E'lp3/2 

^'lpl/2 

Eld5l2 

E 23 I /2 

E-ld3/2 

9 

65.8743670393 

57.3047488634 

56.3996551655 

47.4029580598 

42.0645085975 

45.2829901985 

21 

65.8810348458 

57.3174942083 

56.3981447843 

47.4224876123 

42.1861879711 

45.2877368920 

29 

65.8810254734 

57.3174829474 

56.3980580463 

47.4224848157 

42.1860824080 

45.2876279606 

41 

65.8810257955 

57.3174845038 

56.3980562850 

47.4224891983 

42.1860839907 

45.2876260664 

49 

65.8810258097 

57.3174845532 

56.3980563048 

47.4224893238 

42.1860840610 

45.2876261362 

61 

65.8810258137 

57.3174845640 

56.3980563150 

47.4224893479 

42.1860841046 

45.2876261556 

73 

65.8810258146 

57.3174845663 

56.3980563177 

47.4224893526 

42.1860841130 

45.2876261605 

81 

65.8810258148 

57.3174845667 

56.3980563182 

47.4224893535 

42.1860841146 

45.2876261614 

93 

65.8810258149 

57.3174845669 

56.3980563185 

47.4224893539 

42.1860841154 

45.2876261619 

101 

65.8810258149 

57.3174845670 

56.3980563186 

47.4224893540 

42.1860841156 

45.2876261620 

109 

65.8810258149 

57.3174845670 

56.3980563186 

47.4224893541 

42.1860841157 

45.2876261621 

125 

65.8810258149 

57.3174845670 

56.3980563186 

47.4224893541 

42.1860841158 

45.2876261621 

145 

65.8810258149 

57.3174845670 

56.3980563186 

47.4224893541 

42.1860841158 

45.2876261622 


Table Ic: Same as Table lb but for 4*** order Lagrange FEM. 


5. Order 


T^nod 

Eis1/2 

Elp3l2 

£-1^1/2 

Eld3/2 

E2SI/2 

Eld3/2 

11 

65.8784571165 

57.3086583184 

56.3940420670 

47.4017253841 

42.1771582984 

45.2818654596 

21 

65.8810087294 

57.3174496620 

56.3980829394 

47.4224368621 

42.1858924489 

45.2876670981 

31 

65.8810255741 

57.3174839393 

56.3980557441 

47.4224881501 

42.1860830079 

45.2876254614 

41 

65.8810258056 

57.3174845338 

56.3980562769 

47.4224892693 

42.1860841170 

45.2876260714 

51 

65.8810258150 

57.3174845669 

56.3980563285 

47.4224893538 

42.1860841175 

45.2876261744 

61 

65.8810258146 

57.3174845661 

56.3980563217 

47.4224893521 

42.1860841144 

45.2876261671 

71 

65.8810258149 

57.3174845668 

56.3980563189 

47.4224893537 

42.1860841158 

45.2876261626 

81 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table Id: Same as Table Id but for 5**“ order Lagrange FEM. 
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6. Order 


'^nod 

Eis1/2 

Elp3/2 

^'lpl/2 

Eld5l2 

£'2s1/2 

Eld3/2 

7 

65.8637153942 

57.2828409504 

56.3605696991 

47.3461657413 

41.7817417812 

45.2297330851 

13 

65.8808951903 

57.3169102490 

56.3973587036 

47.4203213172 

42.1844043482 

45.2854747005 

19 

65.8810243462 

57.3174818080 

56.3980718346 

47.4224854807 

42.1860829529 

45.2876435014 

25 

65.8810248977 

57.3174812571 

56.3980529004 

47.4224808348 

42.1860825374 

45.2876204571 

31 

65.8810257970 

57.3174845238 

56.3980563858 

47.4224892554 

42.1860840306 

45.2876261153 

37 

65.8810257899 

57.3174845033 

56.3980564208 

47.4224892313 

42.1860839889 

45.2876263281 

43 

65.8810258112 

57.3174845559 

56.3980563261 

47.4224893294 

42.1860841068 

45.2876261751 

49 

65.8810258148 

57.3174845666 

56.3980563187 

47.4224893531 

42.1860841163 

45.2876261622 

55 

65.8810258149 

57.3174845670 

56.3980563189 

47.4224893542 

42.1860841157 

45.2876261625 

61 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841158 

45.2876261623 

67 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841158 

45.2876261623 

73 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table le: Same as Table Id but for 6**' order Lagrange FEM. 


7. Order 


f^nod 

E'lsl/2 

Elp3/2 

£-1^1/2 

Eld5/2 

£'2s1/2 

E\d3/2 

8 

65.8665180817 

57.2917538980 

56.3821989755 

47.4049510540 

42.0593121642 

45.2786107904 

15 

65.8808831962 

57.3171483037 

56.3981295386 

47.4219539180 

42.1850478466 

45.2877364182 

22 

65.8810247273 

57.3174817834 

56.3980551153 

47.4224839062 

42.1860817524 

45.2876261615 

29 

65.8810255033 

57.3174838576 

56.3980570146 

47.4224881798 

42.1860821285 

45.2876272423 

36 

65.8810258102 

57.3174845542 

56.3980563460 

47.4224893286 

42.1860840744 

45.2876262005 

43 

65.8810258149 

57.3174845669 

56.3980563190 

47.4224893538 

42.1860841147 

45.2876261627 

50 

65.8810258149 

57.3174845670 

56.3980563186 

47.4224893540 

42.1860841156 

45.2876261621 

57 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841158 

45.2876261622 

64 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841158 

45.2876261622 

71 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table If: Same as Table le but for 7^^ order Lagrange FEM. 
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Table Id displays eigenvalues which have been calculated with finite elements of 5*** order. At 76 
mesh points 12 digits have stabilized for all 6 eigenvalues. At 41 mesh points the precision is already 
as good as the precision in Table la at 600 mesh points. A comparison of the eigenvalues in Table 
Id with results of a 6**“ order FEM calculation in Table le shows that a further increase of the order 
leads to a rather weak reduction of the number of required mesh points. At least 73 mesh points 
are necessary in 6*^ order for a precision of 12 digits. As shown in Table If, the reduction of the 
number of mesh points is even weaker when the order is increased from 6*^ order to 7**' order. In the 
subsequent Tables 2a to 2f, results of corresponding calculations with B-spline finite elements are 
shown. In Table 2a, neutron single particle eigenvalues are listed which have been calculated with 
the new B-spline FEM code. A comparison of the numbers with those listed in Table la shows that 
they are identical for equal mesh point numbers. For increasing order of the B-splines, the number of 
required mesh points to obtain a certain precision reduces very similarly to the trend observed in the 
Tables la to If. A comparison of the Tables 2a to 2f with the corresponding Tables la to If shows 
that roughly half the number of mesh points is required in a B-spline FEM in order to achieve the 
precision of a corresponding calculation with Lagrangian finite elements. In Table lb, full precision 
is achieved at 60 mesh points while 121 mesh points were necessary in Table lb. In a calculation 
with 4**' order B-spline elements, 45 mesh points are required as shown in Table 2c whereas 145 
mesh points are necessary with Lagrange elements (Table Ic) to obtain a precision of 12 digits. In 
the 5*^ order B-spline FEM, 34 mesh points have been used (Table 2d) while a corresponding 5*^ 
order Langrange FEM required 76 mesh points (Table Id). The 6*^ order B-spline FEM (see results 
in Table 2e) leads still to a considerable relative reduction of the number of mesh points from 34 to 
30 at the same level of precision while in the 7**' order method still 29 mesh points were required 
(Table 2f). The results shown in the Tables la to If and in the Tables 2a to 2f lead to the conclusion 
that the B-spline FEM has its optimum at 6*^ order whereas the optimal order of the Lagrange 
FEM is at 5*^ order. However, the optimal order may depend on the required precision. 


1. Order 


'^nod 

EisI/2 

-Elp3/2 

Elpl/2 

£■11*5/2 

£'2s1/2 

£■11*3/2 

10 

65.88844974 

57.33594611 

56.42536540 

47.45864481 

42.33620819 

45.34520938 

20 

65.88133752 

57.31827220 

56.39907816 

47.42408909 

42.19281200 

45.28981141 

30 

65.88108502 

57.31763456 

56.39825020 

47.42279518 

42.18737854 

45.28804244 

40 

65.88104432 

57.31753150 

56.39811692 

47.42258515 

42.18648990 

45.28775645 

50 

65.88103335 

57.31750369 

56.39808101 

47.42252842 

42.18624969 

45.28767927 

60 

65.88102944 

57.31749377 

56.39806819 

47.42250815 

42.18616380 

45.28765171 

70 

65.88102777 

57.31748952 

56.39806272 

47.42249948 

42.18612708 

45.28763993 

80 

65.88102696 

57.31748747 

56.39806006 

47.42249529 

42.18610928 

45.28763422 

90 

65.88102653 

57.31748638 

56.39805866 

47.42249306 

42.18609982 

45.28763119 

100 

65.88102628 

57.31748575 

56.39805785 

47.42249178 

42.18609441 

45.28762946 

150 

65.88102591 

57.31748480 

56.39805662 

47.42248983 

42.18608615 

45.28762681 

200 

65.88102584 

57.31748464 

56.39805641 

47.42248951 

42.18608476 

45.28762637 

300 

65.88102582 

57.31748458 

56.39805634 

47.42248938 

42.18608424 

45.28762620 


Table 2a: Same as Table la but for B-spline FEM. 
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3. Order 


T^nod 

£'1s1/2 

Elp3/2 

Elpl/2 

Eld5/2 

E 2 SI /2 

Eld3/2 

10 

65.88136334 

57.31813399 

56.40016227 

47.42348454 

42.18920652 

45.29093938 

20 

65.88102627 

57.31748559 

56.39805814 

47.42249117 

42.18608805 

45.28762926 

30 

65.88102583 

57.31748459 

56.39805636 

47.42248941 

42.18608422 

45.28762623 

40 

65.88102582 

57.31748457 

56.39805632 

47.42248936 

42.18608413 

45.28762617 

50 

65.88102582 

57.31748457 

56.39805632 

47.42248935 

42.18608412 

45.28762616 

60 

65.88102581 

57.31748457 

56.39805632 

47.42248935 

42.18608412 

45.28762616 

70 

65.88102581 

57.31748457 

56.39805632 

47.42248935 

42.18608412 

45.28762616 


Table 2b: Same as Table lb but for B-spline FEM. 


4. Order 


TT'nod 

Eis1/2 

Elp3/2 

£-1^1/2 

Eld5/2 

E 2 SI /2 

E\d3/2 

6 

65.88135693 

57.31850571 

56.40626398 

47.42273069 

42.18460453 

45.29173241 

7 

65.87919142 

57.31238878 

56.39343856 

47.41226119 

42.18527059 

45.28088247 

8 

65.88036297 

57.31616474 

56.40144465 

47.42037951 

42.18049023 

45.29330556 

9 

65.88118037 

57.31768604 

56.39894718 

47.42267251 

42.18838460 

45.28845841 

10 

65.88085949 

57.31705080 

56.39814878 

47.42162032 

42.18564082 

45.28786499 

11 

65.88104483 

57.31749502 

56.39864079 

47.42246069 

42.18618288 

45.28849993 

20 

65.88102585 

57.31748464 

56.39805656 

47.42248949 

42.18608449 

45.28762653 

25 

65.88102582 

57.31748457 

56.39805633 

47.42248937 

42.18608414 

45.28762618 

30 

65.88102582 

57.31748457 

56.39805632 

47.42248936 

42.18608412 

45.28762616 

35 

65.8810258150 

57.3174845672 

56.3980563189 

47.4224893544 

42.1860841164 

45.2876261626 

40 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893542 

42.1860841160 

45.2876261623 

45 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 

50 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table 2c: Same as Table Ic but for B-spline FEM. 


5. Order 


T^nod 

Eis1/2 

Elp3/2 

£-1^1/2 

Eld5/2 

E 2 SI /2 

E\d3/2 

7 

65.88022738 

57.31648520 

56.40099275 

47.42158367 

42.17900675 

45.29112016 

10 

65.88103149 

57.31749002 

56.39852712 

47.42249777 

42.18607033 

45.28816996 

20 

65.88102581 

57.31748455 

56.39805651 

47.42248931 

42.18608408 

45.28762647 

30 

65.8810258149 

57.3174845671 

56.3980563188 

47.4224893542 

42.1860841160 

45.2876261624 

31 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893542 

42.1860841160 

45.2876261623 

32 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893542 

42.1860841159 

45.2876261623 

33 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893542 

42.1860841159 

45.2876261622 

34 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table 2d: Same as Table Id but for B-spline FEM. 


21 

































































6. Order 


T^nod 

Eis1/2 

Elp3/2 

^'lpl/2 

Eld5l2 

E 23 I /2 

E\d3l2 

10 

65.88092534 

57.31722378 

56.39803360 

47.42198939 

42.18580853 

45.28768763 

15 

65.88102458 

57.31748100 

56.39805735 

47.42248186 

42.18608340 

45.28762800 

20 

65.8810258052 

57.3174845391 

56.3980563477 

47.4224892980 

42.1860841341 

45.2876262023 

25 

65.8810258149 

57.3174845669 

56.3980563193 

47.4224893539 

42.1860841164 

45.2876261632 

26 

65.8810258148 

57.3174845668 

56.3980563194 

47.4224893537 

42.1860841155 

45.2876261633 

27 

65.8810258149 

57.3174845671 

56.3980563188 

47.4224893542 

42.1860841159 

45.2876261625 

28 

65.8810258149 

57.3174845670 

56.3980563188 

47.4224893541 

42.1860841159 

45.2876261623 

29 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261623 

30 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table 2e: Same as Table le but for B-spline FEM. 


7. Order 


T^nod 

Eis1/2 

Elp3l2 

^'lpl/2 

Eld5l2 

E 2 SI /2 

E\d3l2 

10 

65.88099742 

57.31745060 

56.39818784 

47.42246110 

42.18582288 

45.28773718 

15 

65.88102496 

57.31748290 

56.39806113 

47.42248693 

42.18607611 

45.28763254 

20 

65.88102580 

57.31748453 

56.39805640 

47.42248929 

42.18608399 

45.28762629 

25 

65.8810258147 

57.3174845665 

56.3980563197 

47.4224893532 

42.1860841144 

45.2876261638 

26 

65.8810258149 

57.3174845670 

56.3980563188 

47.4224893541 

42.1860841159 

45.2876261623 

27 

65.8810258149 

57.3174845669 

56.3980563188 

47.4224893539 

42.1860841157 

45.2876261624 

28 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841158 

45.2876261622 

29 

65.8810258149 

57.3174845670 

56.3980563187 

47.4224893541 

42.1860841159 

45.2876261622 


Table 2f: Same as Table If but for B-spline FEM. 
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log. absolute error log. absolute error 


To complete this study, I repeated the calculations for a large number of mesh points with both 
methods from order to 8 *^ order hnite elements. In the hgures Fig. 8 a to Fig. 8 d, the logarithmic 
errors with respect to the highest precision are shown. Fig. 8 a displays the averaged error taken 
over all 6 single particle energies which have been calculated with the B-spline FEM. In Fig. 8 b, 
these data have been smoothed by taking in addition the average over two neighboring mesh point 
numbers. In the region of precision ranging from 10~^ to both hgures display an enormous 

reduction of the errors for increasing hnite element order up to 5*^ order. Finite elements of higher 
order do not essentially improve the precision in the considered range but entail a higher numerical 
cost. They may improve the precision beyond the error range of [l0~^, . However, precisions 

in the range of [l0~^, are sufficient in most applications. In Fig. 8 b there is an indication 

for 6 *^ order to become optimal order at precisions better than 10“^°. This is in agreement with 
the conclusion that has been drawn form the data in Table 2a to 2f. 

In Fig. 8 c and Fig. 8 d. results that correspond to Fig. 8 a and Fig. 8 b but calculated with the 
Lagrange FEM are displayed. A similar but weaker reduction of the errors is observed for increasing 
hnite element order. A comparison with the results depicted in Fig. 8 a and Fig. 8 b shows that the 
B-spline method requires in general a much smaller number of mesh points than the Lagrange FEM 
in order to provide the desired level of precision. For comparison, in Fig. 8 d, I have inserted those 
graphs of Fig. 8 b which resulted for 5*^ order to 8 *^ order B-spline calculations. 




Fig. 8 c 


Fig. 8 d 




number of mesh points number of mesh points 


Fig. 8: Logarithmic plots of the error which occurs in the B-spline FEM (hgures (a),(b)) and 
in the Lagrange FEM (hgures (c),(d)). The used hnite element orders are indicated. Figs, (a) and 
(c) display average values taken over the six lowest positive neutron single particle eigenvalues. Figs, 
(b) and (d) show the smoothed curves. 
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VII. PERFORMANCE DISCUSSION 


With respect to applications of the above presented B-spline FEM in large scale computations 
with FEM discretizations in more than one dimension, attention should be payed to the performance 
of both methods at the present stage. Therefore, I investigate and compare the run time for both 
methods, the B-spline FEM and the Lagrange EEM. The following discussion is based on data which 
correspond to the performance of the codes on a DEC Alpha 300 MHz. A gnu compiler has been 
used under UNIX to translate the codes. 

The CPU-time depends essentially on the EEM order and the number of used mesh points. The 
times which are displayed in the Figs. 9a and 9b correspond to a single step in which the Dirac 
equation ^ is solved. This procedure encloses essentially the construction and the solution of the 
generalized eigenvalue problem for one Ac-value. It has to be repeated for each k in the solver for 
neutron and proton states and this again over the whole self-consistent iteration. The CPU-time 
which is required for the solution of the meson field equations lies below one percent of that for 
the nucleon states and is therefore neglected. Fig. 9a displays the CPU-time for 5**' order FEM 
as a function of the number of mesh points. Results which have been obtained for other orders 
are almost identical with those shown in the figure. The solid curve displays CPU-times resulting 
from Lagrange type finite element discretizations while the dashed curve has been obtained with the 
B-spline EEM. An explanation for the higher numerical cost in applications of the B-spline method 
is given by Figs. 4 showing that more matrix elements have to be calculated in the case of B-spline 
FEM. 

However, as demonstrated in the Figs. 8, the number of required mesh points for equal numerical 
precision is half of that in the Lagrange FEM. At equal numerical precision, the solid curve has 
to be compared with the dot-dashed curve of a B-spline calculation. This clearly demonstrates an 
enormous reduction of the numerical cost for the B-spline FEM. 

In Fig. 9b, CPU-times are plotted as a function of the FEM order and compared for both methods. 

The precision of the numerical solution which has been obtained with the B-spline FEM (dashed 
line) is much higher that that obtained with the Lagrange FEM (solid line). At equal numerical 
precision, one should compare values of the dashed line with values of the solid line at double order. 

Fig. 9a Fig. 9b 




Fig. 9: CPU time as a function of the number of used mesh points (a) and as a function of 
used order (b) for both methods, B-spline FEM and Lagrange FEM. 


VIII. PROGRAM STRUCTURE 

The program is coded in C-F-F. The implementation of the relativistic mean field model in the 
Hartree approximation for spherical doubly-closed shell nuclei has been described in Ref. |^. In 
this section we only describe the changes that have been made in order to modify the program to 
B-spline techniques. 

The main part of the program consists of seven classes: MathPar: numerical parameters used 
in the code, PhysPar: physical parameters (masses, coupling constants, etc.), FinEl: finite ele- 
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ments, Mesh: mesh in coordinate space, Nucleon: neutrons and protons in the nuclear system, 
Meson: mesons and photon with corresponding mean fields and the Coulomb field, and the class 
LinBCGOp. A detailed description of these classes can be found in Ref. Q]. 

Two new classes FinElBsp and BSpline have been added to the code. The implementation 
can be found in the source files finelbsp.cc, bspline.cc and the corresponding header files finelbsp.h, 
bspline.h. The class FinElBsp contains the following new member functions: 

FinElBsp 0; 

FinElBspO; 

void EinElBsp::alloc ( int ord ); 

void EinElBsp::free(); 

void EinElBsp::make( int ord ); 

double FinElBsp::n( int Hoc, int ife, int iga, int I, bool zero ) const; 

double FinElBsp::dn( int Hoc, int ife, int iga, int I, bool zero ) const; 

double FinElBsp::func( double const* u, int ife, int iga ) const; 

double FinElBsp::func( double const* u,int ife,int iga,int I,bool zero ) const; 

void EinElBsp::eval(); 

inline int FinElBsp::order() const; 

inline int FinElBsp::nloc() const; 

inline double FinElBsp::n( int Hoc, int iga ) const; 

inline double FinElBsp::dn( int Hoc, int iga ) const; 

inline double FinElBsp::operator()(double const* u,int ife,int iga) const; 

inline double FinElBsp::operator()(double const* u,int ife,int iga,int I, bool zero ) const; The two 
methods EinElBsp() and FinElBsp() describe the constructor and destructor of objects (B-spline 
finite elements) of the class. The method make( int ord ) provids the B-spline reference element 
with data (amplitudes) on the shape functions and their derivatives. Access to these data is given 
through the methodes n( int Hoc, int iga ) and dn( int Hoc, int iga ^ . In a first step, make( int ord ) 
allocates memory for the shape functions using alloc( int ord ). In a second step, the method eval() 
is called generating the amplitudes of the shape functions through an operator of class BSpline 
. The overload member functions n( int Hoc, int ife, int iga, int I, bool zero ) const and dn( int 
iloc, int ife, int iga, int I, bool zero ) const are used in the calculation of the stiffness matrices. 
They take into account boundary conditions. The method func( double const* u, int ife, int iga ) 
const provides the interpolated amplitude of solutions on the Gauss submesh in any finite element of 
global index ife. The overloaded version/iinc(^ double const* u, int ife, int iga ) const takes boundary 
conditions into account. 

In the class BSpline, the following methods are implemented BSpline( int ord ) 

B Spline 0 

operator()( double const* p, double& f, double& df, double x ) 
inline int BSpline::order() const; 

BSpline( int ord ) and BSpline() describe the constructor and the destructor of the class. An 
operator is used to carry out the B-spline algorithm whenever access to B-spline amplitudes is 
requested through a call of an object of the class with corresponding arguments. 

The organization in the construction of stiffness matrices in other parts of the code has been 
changed accordingly. However, the essential structure has been maintained so that quick changes 
for applications of Langrange type finite elements (defined in class FinEl) are possible. An essential 
difference roots in the relation between number of nodes on the global mesh and the number of finite 
elements which is given by n"°‘^ = -|- In the version using Langrange type elements this 

relation is n“°‘^ = -|- 1. 

For the diagonalization of the generalized eigenvalue problems the bisection method has been 
replaced by a combined Cholesky decomposition and householder method which is slower but allows 
for a higher precision. A heapsort algorithm orders eigenvalues and eigenvectors. The routines and 
the eigensolver are implemented in the source file eigen.cc and the header file eigen.h. 
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APPENDIX: APPENDIX 


For the FEM discretization of the Klein-Gordon equations we use the ansatz 


<^{p) = ^(PpBp{r) 

(Al) 

E 

O 

(A2) 

= '^P°pBp{r) 

P 

(A3) 

A°{r) = Y.AlBp{r). 

(A4) 


p 


where the node variables CTp, cUp, Pp and A® correspond to field amplitudes at the mesh point p. For 
the Klein-Gordon equations we use the same type of shape functions Bp(r), and the same mesh as 
in the FEM discretization of the Dirac equation (p^). Using again the method of weighted residuals 
with test functions Wp{r) = r^Bp{r), the following algebraic equations are obtained 


P 

^(u)p/(r) 

P 

P 

^(wp/ (r) 


o2 2 /(Z -b 1) 2 

-dr - dr -i-5-h 


— dr - dr + 

r r 


1(1 + 1) 


+ 


p,2 2 Z(Z -b 1) 2 

— dr — -dr H- ^ -b rup 


Sa($i(r),...,$A(r)) 


Sa;($i(r),...,$A(r)) 


Sp{^i{r),...,^A{r)) 


-d^r--dr+^-H^ 


Bp{r)'jap = (wp,{r) 
Bp(r)^tjp = (wp,{r) 
Bp{r)'jpl = (wp>{r) 


Bp{r)'^A°p = (wp,{r) sc('hi(r), 


The resulting matrix equations read 


(A5) 

(A6) 

(A7) 

(A8) 


5*1 + /(/ + 1 ) • 5*2 + • 5*3 + 5*4 

5*^ +/(/+ 1) • 52^ + 5^ 


- As) 

(j = ' 


-♦0 -^v) 

■ UJ = 7^ ^ 


S? + l{l + l)-S^ + ml-S^ 


aO , ^ A 0 

sf +l{l + l)-S^ 


-0 -A3) 

■ p = ' 


■A' 


'0 _ Aem) 


The node variables o-p, Pp and Ap are grouped into the vectors a = (cri, cu^ = 
p° = A° = (A?,...,A°)^, and 


9^ -b 2 r ^ dr 


Sr = Sr = = Sf = (^Wp,{r) 

= SP = si = (^Wp,ir)\rABp{r) 
Si = Si = Si = si = ^«;y(r)|Bp(r)^, 


Bp{r) 


Si = (^Wpi (r) 


520-(r) -b530-(r) 


Bp{r)). 


The components of the right hand side vectors are defined as 

_ 


= -g^(^Wp,ir) ps(r)^, 
= 9<^(wp^{r) pv(r)^, 

rfi =gp(wp>{r) P3(r)^, 


PpA = e(wp,{r) pem(r)y 


(A9) 

(AlO) 

(All) 

(A12) 


(A13) 

(A14) 

(A15) 

(A16) 


(A17) 
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The nonlinear equation for the (T-field is solved in an iterative procedure. The nonlinear terms are 
included in the global stiffness matrix (see S 4 in Eq. (A9)). In the iterative solution matrix elements 
that contain nonlinear terms are calculated using the field a{r) obtained in the previous iteration 
step. 
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Fig. 3 e 
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